Density profiles and voids in modified gravity models 
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We study the formation of voids in a modified gravity model in which gravity is generically stronger 
or weaker on large scales. We show that void abundances provide complementary information to halo 
abundances: if normalized to the CMB, models with weaker large-scale gravity have smaller large 
scale power, fewer massive halos and fewer large voids, although the scalings are not completely 
degenerate with erg. Our results suggest that, in addition to their abundances, halo and void 
density profiles may also provide interesting constraints on such models: stronger large scale gravity 
produces more concentrated halos, and thinner void walls. This potentially affects the scaling 
relations commonly assumed to translate cluster observables to halo masses, potentially making 
these too, useful probes of gravity. 



I. INTRODUCTION 

A wealth of observations, from WMAP, supernovae la, 
galaxy clustering on large scales, and cross-correlations 
between galaxies and the CMB, suggest that we live in a 
spatially flat Friedmann-Robertson- Walker universe cur- 
rently dominated by either a cosmological constant or re- 
pulsive dark energy. The best fit for the dimcnsionless en- 
ergy density parameters are fi m = 0.28 and Q\ = 1— f2 m , 
within the concordance ACDM Model Tj. However, be- 
cause this requires the vast majority of the energy density 
in the universe to be in two unknown substances, dark 
matter and dark energy, there is considerable interest in 
alternative interpretations of what the data imply. 

The key equation of General Relativity, Einstein's 
equation, relates the curvature and the expansion rate 
of the Universe to its matter and energy content. The 
current paradigm is to modify the matter content of the 
universe, by including dark matter and dark energy, to 
account for observations. Instead, however, we might 
modify how the universe curves in response to mat- 
ter, which would mean modifying our theory of grav- 
ity. There are many wa ys t hat one could modify gr avity 
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but the one that we will focus on in this paper is purely 
phenomenological. The assumption is that it is merely 
a modification to the potential, changing the form from 
that of standard gravity to: 



(r) = Gm 



1 + a(l - e- r ' r -) - a(l - e~ r / r ") 



(1) 



Thispotential (with r c — > oo has been studied before [13, 
t^HHHHLIl^l and in some sense, is an interesting limiting 
case of some f(R) models 0, [25j]. In what follows, we 
will use r c r s but finite, following [23[. The precise 
choice of r c does not matter for any of our conclusions. 

In [23[ we studied the formation of nonlinear objects 
in this model - the first study of virialized dark mat- 



ter halo abundances in any modified gravity model. In 
Section [Til we study the evolution of the structures that 
are in some ways the opposites of halos, namely voids. 
Along the way, we show that the density profiles of ob- 
jects in these models exhibit interesting departures from 
standard gravity models. Section IIIII shows how we use 
insights from the evolution model to estimate void abun- 
dances. Section [TV] shows that our model captures the 
essence of how void abundances depend on a and r s in 
numerical simulations of structure formation with this 
modified potential. A final section summarizes our re- 
sults. An Appendix shows the corresponding trends for 
halo profiles; these potentially allow X-ray observations 
to provide powerful probes of modifications to gravity. 



II. EVOLUTION OF UNDERDENSITIES 

This section describes the evolution of initially un- 
derdense spherically symmetric regions for the modified 
gravity model described in the previous section. 

A. The spherical tophat model 

Following flfj, we consider the evolution of a spheri- 
cal patch of density different from the background in an 
expanding universe. The spherical collapse calculation 
begins with the statement that the force driving the ac- 
celeration is related to the gravitational potential by 



dt 2 



= F = -V$(r). 



This can be integrated once to get 



•*(r)=<7, 



(2) 



(3) 



where $(r) is the integral of the potential over the mass 
distribution, and C is the total energy of the patch, which 
is constant. In standard gravity, the potential of a shell 
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center of the sphere, so $(r) reduces to GM(< r)/r. 
The constant C can be related to the initial overdensity 
and/or expansion rate of the patch: the initial expan- 
sion rate is given by the Hubble expansion rate of the 
background in which the patch is embedded, namely in 
comoving coordinates ii — 0, so fi — 04X1 — hi / aiifliXi) , 
so {dr/di)i = HiVi. Including a cosmological constant 
presents no conceptual difference. 

In standard gravity one can directly solve this equa- 
tion. The solution is a cycloid, which is usually writ- 
ten in p arametric form involving sines and cosines (but 
see (27| for simple but accurate approximation). Only 
objects which are sufficiently overdense initially will col- 
lapse (meaning the size of the patch has become vanish- 
ingly small) at the present time. The critical value of the 
initial overdensity required for collapse today , <L , docs 
not depend on the initial size of the patch [26, 28]. This 
scale independence of S c is a result of Birkhoff 's Theo- 
rem: the evolution of a tophat sphere is the same as that 
given by Friedmann's equations, so the actual size of the 
patch drops out. (Indeed, for triaxial perturbations, S c 
depends on the size and shape of the initial patch [291].) 

When gravity is modified, things are no longer so sim- 
ple. For example, when the potential is given by equa- 
tion (fT]), Birkhoff's theorem no longer applies: A particle 
offcenter in a uniform spherical shell will feel a force from 
the shell because the force no longer varies as 1/r 2 . This 
has two consequences. First, equation © can still be 
integrated once to get equation ©, only now <i>(r) has 
contributions from both the internal and external mass 
distributions. We can get $ by integrating equation (TT|) 
of the mass distribution, about which more shortly, leav- 
ing C and the initial value for dr/dt to be determined. 
As before, C is the total energy (constant in time), and 
we set (dr/dt)i = HiTi. 

Second, whereas evenly spaced concentric shells remain 
evenly spaced in the standard tophat model, this is no 
longer the case when the potential is modified. This is 
most easily seen by writing as the sum of three terms: 
that due to the Newtonian 1/r part, and those due to 
the Yukawa term from shells internal and external to the 
point (or shell) of interest (see Appendix A for details). 
If 6 denotes the overdensity of the source shell, then the 
internal and external contributions are both proportional 
to a5, but they have opposite signs. As a result, the 
initial tophat perturbation develops a nontrivial density 
profile (see Appendix [A] for examples) , meaning interior 
shells may cross one another as the initial 'boundary' 
around them collapses. When a > 0, initially overdense 
perturbations become more centrally concentrated than 
when a = 0, as the contribution from the modified part 
of the potential pulls mass away from the boundary in 
both directions (provided the boundary is comparable to 
r s ). When a < 0, then the perturbation develops a ridge 
at its boundary, as mass is pulled towards the boundary 
from both directions. 




Comoving Radius (in Mpc/h) 




Comoving Radius (in Mpc/h) 




Comoving Radius (in Mpc/h) 

FIG. 1: Evolution of the density profile of a void. The three 
lines are the density at the initial time, the halfway time (what 
would correspond to the turnaround time in halo formation) , 
and the final time. Of note is that ridge formation occurs in 
all three cases, though for negative a the density is noticeably 
enhanced outside even the ridge. 
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B. Density profiles 

For voids, much of the preceding discussion remains 
valid. In standard gravity (a = 0), the solution is not a 
cycloid, but similar: the sines and cosines are replaced 
by their hyperbolic counterparts (however, the simple ap- 
proximation of [27j remains accurate even in this regime) . 
When a ^ 0, then the previous discussion continues to 
apply. However, because the sign of 5 is reversed, the de- 
pendence on a is also reversed. As a result, mass moves 
towards the boundary of the profile from both sides when 
a > 0, but away from it when a < 0. 

Figure [T] shows the evolution of a number of concen- 
tric shells within which the initial density was uniform, 
but slightly below that of the background. We used the 
methods of [23[ to calculate the evolution of the shells. 
(Note, however, that whereas few shells are needed to ac- 
curately estimate the evolution of the boundary of the ini- 
tial tophat perturbation, many more are needed to study 
the evolution of the density profile.) When a = 0, then 
the expansion of the inner shells pushing against the ex- 
terior shells leads to the formation of a ridge [3(| and 
references therein]. This remains true when a ^ 0, al- 
though, when a > then the ridge is smaller in physical 
size, whereas the opposite is true for negative a. In par- 
ticular, for the a < case, the ridge extends further and 
trails off much more slowly. In addition, when a / 0, 
then the density in the interior region does not remain 
constant - this is clearly more evident in the negative a 
case, but also evident in the positive a case. Whereas 
the depedence of the interior density profile on a may be 
difficult to detect in observations, the effect on the ridge 
may be observable. Of course, to really exploit this ef- 
fect, one must study more realistic initial density profiles 
- this is beyond the scope of our work. 



C. Critical underdensity for void formation 

One sense in which underdensities are more difficult to 
study than overdensities is that the concept of a criti- 
cal time, the analogue of the collapse time, is not well- 
defined. While it is clear that a halo is collapsed at 
that time when the radius of its outermost shell reaches 
zero, what is the appropriate condition for undersities? 
For standard gravity, one uses the condition of "shell- 
crossing" which occurs when initially interior shells first 
cross initially exterior shells: the objects for which is true 
are about 0.2 times denser than the background [3(J. Be- 
cause of the ambiguity associated with shell crossing in 
modified gravity models, we will use this value, an un- 
derdensity of —0.8, to define voids even when a =/= 0. 
This critical nonlinear underdensity corresponds to some 
critical initial underdensity, 5 V . When a = 0, S v ph 2.86 
independent of the size of the underdensity, for the same 
reason that 8 C w 1.686 is independent of mass. However, 
when a ^ 0, we know S c depends on halo mass (23|, and 
so we expect S v to depend on void volume. 
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FIG. 2: Ratio of initial density required for void (solid) and 
halo (dashed) formation at the present time to that in the 
standard gravity model, when the background cosmology is 
ACDM, r s = 5/i _1 Mpc and r c = 70/i _1 Mpc. From top to 
bottom, curves show models in which a — —1, —0.5,0.5 and 
1 (note that a — is standard gravity). 



To find this dependence, we study the evolution of 
underdense patches starting from a large grid of initial 
underdensities and sizes. Once we have the scale de- 
pendence of S v (Ri), we can relate it to a mass scale 
by M = (4ir/3)Rfpi(l + Si) « (47r/3)i?f, the last rela- 
tion holding because Si is small. Thus, we can calculate 
8 V {M). Intuitively, we expect that if, in the course of its 
evolution, the size of a patch never exceeds r s , then we 
expect the large-scale modification to gravity is inconse- 
quential. Hence, for sufficiently small voids, which form 
from smaller patches in the initial field, we expect S v to 
be the same as in standard gravity. For large voids, the 
effect of the modification should be stronger. For posi- 
tive a voids should be easier to form, so we expect \S V \ 
to be smaller, whereas for negative a the opposite should 
be true. 

The solid curves in Figure show that <5 V (M) has the 
expected dependence on M and a. Note that we have 
chosen to express the scale dependence of S v in terms of 
M, for ease of comparison with the scale dependence of 5 C 
shown in Fig. 2 of [23| , and reproduced here as the dashed 
curves. Of course, we could have expressed it in terms of 
the initial size Ri , or in terms of the final size r v . This 
is because we have defined voids as being 0.2 times the 
density of the background, making r v — 5 1 / 3 Ri ~ 1.7 R4: 
the comoving radius of a void is 1.7 times larger than it 
was initially. A close inspection of the figure will show 
that S v departs from its a = value at slightly smaller 
mass scales than does S c . This is because voids expand, 
so smaller mass scales can eventually exceed r s in size, 
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and so notice that a ^ 0. 

III. STATISTICS OF VOID ABUNDANCES 

When combined with the spherical evolution model, 
excursion set methods p8l. l3ll I32L l33l. l34j are commonly 
used to study the abundance of virialized objects. They 
can easily be extended to the study of voids [30] ■ In short, 
this is done by relating the properties of such objects with 
the initial (rather than linearly evolved) density fluctua- 
tion field [23| . 



A. The excursion set method 

Now, as shown, when the potential is modified, then 
5 V is no longer scale- independent. Because it depends 
on mass, the relevant excursion set problem is one with 
a 'moving' rather than 'constant' barrier, so it is of the 
type first studied by [33| . Figure [3] shows the result of 
using this formalism to estimate the abundance of voids. 
Briefly, making this estimate requires that one generate 
an ensemble of random walks in the (Si, Si) plane, where 
Si = u|(M) is the variance in the initial fluctuation field 
when smoothed on scale r$. Since (Ji(M) is a monotonic 
function of M , the variables Si, M and r 4 are essentially 
equivalent to one another. In particular, 

s ^jdkk^ w2{kril (4) 

where W(x) — (3/x 3 ) (sinx — xcosx). For halos, one 
then finds the first crossing distribution f(Si)dSi of the 
'barrier' S C i(M) = S C i(Si). The abundance of objects is 
given by 

^M-Lfvw*. (5) 

For voids all this is slightly more complicated. The 
relationship between f(Si) and dn/dhiM dlnM = 
(p/M) f(Si)dSi is unchanged, but one must account for a 
subtlety in the agrument. To estimate halo abundances, 
one is careful to choose the first upcrossing of the collapse 
barrier, thus giving the largest smoothing radius which 
is sufficiently overdense that it will collapse today. The 
larger scale environment is irrelevant: for example, if the 
larger scale density is negative, the end result would just 
be a halo which formed inside a void. With voids, though, 
one must worry about underdense regions which are sur- 
rounded by large scale overdensities which would have 
collapsed by the present time - the void-in-cloud prob- 
lem. Namely, if on some scale the random walk crosses 
below the critical density for void formation, but on a 
larger scale it crossed above the barrier for halo forma- 
tion, then this would be a void which was crushed as the 
region around it collapsed. Therefore, such walks should 
not contribute to the estimated void abundance. As a re- 
sult, rather than a one barrier problem, we instead have 




lg V (in (Mpc/h) 3 ) 

FIG. 3: Mass function of voids. Black is for standard gravity 
(a — 0), red is for a = —1, blue for a = 0.5. In this case, all 
three mass functions are calculated assuming that the initial 
power spectrum is the same, meaning that the power spec- 
trum at redshift zero is different. 



a two barrier problem [3(ij |: S v must be crossed before S c , 
is crossed. (A more conservative version of this argument 
sets the barrier associated with the collapsing object to 
equal that for turnaround rather than full collapse.) 

Because there are two barriers, the first crossing distri- 
bution has a peaked shape, with the distribution being 
cut off at both high and low masses, on the high end due 
to the rarity of such extreme underdensities, and at the 
low end, because S c < \6 V \. In large part, we will not 
be interested in the intricacies of which halo formation 
barrier we should use, as this largely affects only the low 
mass end of the mass function, where the modification to 
gravity is smaller. 

Figure [3] shows the resulting void mass functions for 
standard and modified gravity. (We do not show the 
extremely low mass regime, which is most sensitive to 
the void-in-cloud cutoff.) Notice that when a is positive 
there are more large voids and fewer small voids, whereas 
the opposite is true when a is negative. This is from two 
effects: one is that any given walk tends to cross later 
when gravity is weaker (a < 0) as the barrier is lower 
(recall <5„ is negative) and so voids tend to be shifted 
towards lower mass. The second effect is that because it 
is harder to form high mass halos in weaker gravity, some 
walks that cross 5 C in stronger gravity may not cross 5 C 
for weaker gravity, however due to the fact that they are 
near S c at high mass, it is still difficult for such a walk to 
cross 6 V even at low mass. The first effect is likely more 
significant, but both contribute to the difference at low 
mass. 
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B. Initial conditions 

One thing to note about the histograms in figure [3] 
is that the difference between the standard gravity and 
modified gravity void mass functions comes from two 
sources. Perhaps the more obvious is that modifying 
gravity modifies the behavior of individual structures, 
changing how quickly they form. This information is en- 
capsulated in the barrier shape, and hence is accounted 
for using excursion set methods. The result of a shifted 
barrier is a shift in the function f(S), the first crossing 
distribution. The other source is the modification to the 
power spectrum. In the study above, we assumed that 
the initial power spectra were the same for all the models 
- what is usually called CMB-normalization. However, 
because gravity is modified, the late time power spectra 
(in both linear and nonlinear theory) depend on a (and 
r s ). So now there are two other possible cases. We could 
consider modified gravity but with the initial conditions 
chosen so that the z — power spectrum is that of stan- 
dard gravity, or we could consider standard gravity but 
with the initial conditions modified so that the z — 
power spectrum is that of modified gravity. In excursion 
set terms, either of these choices affects the relationship 
between M and S. 

For halos, whether one chooses to match the initial 
power spectrum or the final power spectrum matters [23| , 
because it changes the mapping between S and M (or 
i?), and this mapping affects two steps of the calculation: 
the transformation of S(M) into S(S), and the relation 
between f(S) and n(M)dM. The same is true for voids 
so, in what follows, we will consider some of these other 
possiblities as well. 



IV. COMPARISON WITH SIMULATIONS 

We now compare our spherical evolution predictions 
for void abundances with measurements in the simula- 
tions of [2l[ . These simulations followed the evolution of 
128 3 particles in a periodic box of size 100/i -1 Mpc, for 
various choices of a and r s . In all cases, the background 
cosmology was flat ACDM with fl = 0.3, and the par- 
ticle mass was 1.1 x 10 10 M Q . The a — simulation, 
with standard initial conditions has as = 1.0 at z = 0. 
The corresponding runs for a = 0.5 and a = — 1 have 
<Tg = 1.10 and 0.84 respectively. Following our discussion 
of how the halo and void counts depend on the shape 
and normalization of the initial power spectrum, we also 
study results from a — simulations in which the initial 
power spectrum was modified so that, at z = 0, it has 
the same shape as the two a / cases (erg = 1.10 and 
0.84). 

The simulations were analyzed using the void finder of 
[HI . Figure [4] shows the halos and voids found in three 
runs, all with the same iitial conditions. Notice that the 
voids are largest when a — 0.5 and smallest when a = 
— 1, in qualitative agreement with expectations. Since 




20 40 60 80 100 
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FIG. 4: Spatial distribution of halos and voids in three real- 
izations of the same initial conditions evolved to the present 
time using different values of a. Each panel shows voids and 
halos whose centers lie in a slice that is 6/i -1 Mpc thick: top 
panel shows halos and voids in a standard gravity run (a — 0) , 
middle is for a = 0.5, and bottom for a — —1.0. The halos 
become more strongly clustered and the voids larger as a in- 
creases. 
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FIG. 5: Ratio of mass function of voids to that in standard 
gravity, when the initial power spectrum is the same. Red 
is for a = —1, blue for a = 0.5. The histograms (TO BE 
ADDED) are the results from the simulations, the solid lines 
are the results from theory. 



the a = 0.5 run has more large scale power at z = this 
is not surprising. Figure [5] shows the results of a more 
quantitative comparison: although qualitatively similar, 
the predicted effects (solid lines) for large voids are larger 
than we find in the simulations, and the effect on small 
voids is smaller than is observed. 

We explore the dependence on choice of initial con- 
ditions in two steps. The top panel in Figure [6] shows 
the ratio of void counts in two standard gravity runs, 
one with initial conditions modified to produce the same 
2 = power spectrum as q ^ 0, and the other with the 
standard a% = 1 initial conditions (note that the modi- 
fication to the initial conditions is different for the two 
values of a shown). This plot is qualitatively similar 
to the previous one, because stronger large-scale grav- 
ity (a > 0) is qualitatively like having more large scale 
power, but the differences between these two plots shows 
that the result of modifying gravity is not degenerate 
with modifying the shape and normalization of the ini- 
tial power spectrum. The bottom panel shows this infor- 
mation slightly differently: here, the void counts in the 
a ^ run are ratioed to the counts in the a = run in 
which the initial conditions were tuned to produce the 
same z — power. The fact that the ratios are not unity 
implies that there is more information in the void size 
distribution than in the power spectrum itself. 

Finally, Figure [7] shows the ratio of void counts when 
it is the initial conditions in the the modified gravity 
runs which have been tuned to produce the same z = 
power spectrum (rather than tuning the a — initial 
conditions). Unfortunately, we do not have simulations 
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FIG. 6: Ratio of mass function of voids in standard grav- 
ity with modified ICs to standard gravity with standard ICs. 
Top: The ICs are chosen so that at z = 0, the simulations had 
the same power spectrum as the modified gravity simulations. 
Bottom: ICs for the a — runs have been tuned to produce 
the same power spectrum at z — (this tuning is different for 
the two values of a shown). In both panels, red is for a = — 1, 
blue for a = 0.5. 



of this case, but we again see that the ratio is predicted to 
differ from unity, indicating that modifications to gravity 
are not degenerate with changes to the initial conditions. 



V. CONCLUSIONS 

We studied a particular modification to large-scale 
gravity, which has two free parameters: a scale r s , and 



7 




1 2 3 4 5 

lg V (in (Mpc/h) 3 ) 



FIG. 7: Ratio of void counts in modified gravity simulations 
to those in standard gravity. Here, the standard gravity runs 
use standard initial conditions, but the modified gravity runs 
use different initial conditions, chosen so that the z = power 
spectrum is the same as it would be for standard gravity. Red 
is a — —I, blue is a — 0.5. 



an amplitude a (equation [T]) . Previous work has shown 
that halo abundances can be modeled well enough [23| 
to use them to constrain such models. We present an ex- 
ample of this in Appendix [B] We argued that since voids 
are so large, and almost fill space, they too can be useful 
probes. In particular, the study of the formation of large 
voids provides a useful counterpart to the the study of 
the formation of massive halos because both probe the 
nature of gravity on large scales. This is timely, given 
that large, statistically representative void samples are 
only just becoming available [e.g. l3cl \37\. 

Whether stronger large-scale gravity produces more or 
fewer large voids depends on how one normalizes the 
models. When normalized to have the same power spec- 
trum at early times (the most common choice), stronger 
gravity produces more large voids (Figure [3]). This case 
also produces more massive halos, so the net result is 
qualitatively like having standard gravity with more large 
scale power (Figure |4j. Although there are quantita- 
tive differences which allow one to distinguish between a 
larger normalization and modified gravity (Figures EH?]), 
we note that this fact alone - a mismatch between the 
CMB- and later-time determinations of the amplitude 
of the power spectrum on 10 Mpc scales when standard 
gravity is assumed - are generic signatures of modifica- 
tions to large-scale gravity. Note that this particular sig- 
nature has the same sign for halos and for voids - stronger 
large scale gravity means more massive voids and more 
large voids. In contrast, if gravity is standard but the 
initial conditions were non-Gaussian, then halo and void 



abundances are modified in qualitatively different ways, 
at least for the local non-Gaussian initial conditions of 
current interest [38j . 

However, if normalized to have the same power at 
z = 0, the dependence of void (and halo) abundances 
on whether large-scale gravity is stronger or weaker de- 
pends on how one chooses to do this. If this is done by in- 
creasing (decreasing) the initial power spectrum in mod- 
els with a < (a > 0), then one still expects the stronger 
gravity models to produce larger voids (Figure [7]), al- 
though the predicted abundances are quantitatively dif- 
ferent. However, if one modifies the standard gravity 
initial conditions to match those of the a ^ mod- 
els (increase/decrease initial large scale power to match 
a > 0/a < 0), then one predicts more (fewer) large voids 
when a < (a > 0) relative to the a — case. This 
dependence on how the models were normalized is quali- 
tatively similar to the trends seen for massive halos [23| . 

We presented a method for estimating these effects 
which is in good qualitative greement with the simula- 
tions, but there are quantitative differences. However, 
the agreement is not as good as it was for describing ha- 
los. Larger simulations are required to determine if this 
is due to the relatively small size of the simulation boxes 
available to us, or to some more fundamental problem 
with our analysis (see [2!| [34| for discussion of the draw- 
backs of the excursion set approach). 

We also showed that the density profiles of voids and 
halos may also provide interesting probes of modified 
gravity. When the initial profile is a tophat, then the 
evolved halo profile in the a < case generally has a cusp 
at the virial radius (mass flows towards the boundary), 
whereas when a > the halos are more centrally concen- 
trated (Figure [5]). The structure of the void walls also de- 
pends slightly on a (Figure[T|). We provided tentative ev- 
idence that the profiles of halos in simulations do depend 
on a (Figure 0, and hope that these initial measure- 
ment motivate further study of this interesting problem. 
It will be interesting indeed if these trends persist in sim- 
ulations with better mass and force resolution, because it 
is conceivable that effects of this magnitude will soon be 
measured by weak lensing surveys. The more dramatic 
effects associated with initially tophat profiles potentially 
provide an interesting X-ray signature of modified gravity 
- it would be interesting to simulate such models using 
SPH codes. 

In conclusion, we expect our results will prove useful in 
studies which use the large scale distribution of galaxies, 
and the structure of galaxy clusters, to constrain large 
scale modifications of gravity. 
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APPENDIX A: HALO DENSITY PROFILES 

A complementary study to the one in the main text 
regarding the evolution of the density profile of voids is 
a study of the evolution of halo density profiles. Figure[5] 
shows one such study. In this case, the three objects all 
started as tophats with the same initial overdensity, but 
different values of a. The solid and dashed curves show 
the density profile for two different values of r s and at 
three different times, one corresponding to turnaround, 
one to when one shell reaches a density of 200 times the 
background, and the final to when the density enclosed in 
the outer shell reaches 200 times the background. Since 
gravity has a different strength in the three cases, and 
the patches all started with the same overdensity, the 
final collapse happeds at different times. For positive 
a, objects collapse faster, whereas for negative objects 
collapse slower, with the difference decreasing as r s in- 
creases. This is why the effect is noticeable earlier for 
negative a: the turnaround time is later than for stan- 
dard gravity, so the deviation is greater. 

The evolved profiles in the two modified gravity cases 
are strikingly different from standard gravity, with the 
formation of sharp spikes at the edge of the halo in the 
negative a case, and the formation of a large central peak 
in the positive a case. Whether these profiles would hold 
through virialization is unknown, but since the effect is 
so strong during the period before collapse it seems likely 
that the profiles of virialized haloes would be affected. If 
so, then halo profiles might provide new interesting con- 
straints on a. However, to place such constraints, one 
would have to extend our analysis to more realistic ini- 
tial profiles which do not have sharp edges - although it 
is likely that the formation of the cusp at the boundary 
when a < will remain (see Figure [5] and related discus- 
sion). Note that for larger r s , the difference from stan- 
dard gravity is smaller, because the perturbation spends 
essentially all its time at scales which are smaller than 

These trends are perhaps most easily understood by 
writing the potential due to an extended source at scales 
that are beyond its boundary. The potential at r due to 
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FIG. 8: Evolution of the enclosed density of initially tophat 
overdense regions in standard (top) and modified gravity 
(middle and bottom panels show results for a = 0.5 and 
a = —1.0 respectively). In the modified gravity plots, solid 
and dashed curves are for r s — 5/i _1 Mpc and r s — 10/i _1 Mpc 
respectively. In all cases, the mass of the final halo is 
10 15 ' 3 /i -1 Mq, and the density profile is shown (from bottom 
to top) at turnaround, when a shell reaches a local density 
of 200 pb, and when shell within which the enclosed mass is 
10 15 ' 3 /i -1 M© reaches an enclosed density of 200 pb- (For stan- 
dard gravity these last two conditions occur simultaneously.) 
Increasing r s leads to halo profiles that are more similar to 
standard gravity. 




1 2 3 4 5 




1 2 3 4 5 




1 2 3 4 5 



R / R «lr 

FIG. 9: Profiles of the enclosed density for the most massive 
halos in the simulations started from identical initial condi- 
tions but with a = (top), alpha — 0.5 (middle) and a = — 1 
(bottom). 



10 



a top perturbation which contains mass M within radius 
R can be written as the sum of the standard gravity piece 
plus a term which arises from the modification: 



$(r) 



GM GSM r 
a — 



e~ r ' r ' T{R/r.) (Al) 
- r l r °F{R/r c ) 



where <5M = M - p4nR 3 /3 and 



a; cosh(a;) — sinh(x) 



(A2) 



Note that T > 1 for all x. However, T — ► 1 when x -C 1, 
so it is easy to see that when r s ^ R then $ looks just 
like it does standard gravity. If the overdensity profile is 
a power-law of slope —1 within R, so = 2irR 3 , then 
the expression for <!> remains true, but T(x) — > [cosh(a;) — 
l]/(x 2 /2): again, J 7 > 1 for all x, and — > 1 when 
i«l. In practice, we are most interested in the regime 
in which R is smaller than a few times 2r s , for which 
the difference between these form factors is small. This 
means that one can approximate the force outside the 
perturbation by assuming it is a tophat, and not worrying 
about its internal structure. In fact, this suggests that 
one can solve for the evolution of the boundary of the 
perturbation without worrying about the fact that the 
internal structure is changing; in essence this is why the 
analysis in [23[ was so straightforward. 

The evolution of the profile within the perturbation 
is more easily studied using the force, which is V r <&(r), 
and can be split into a Newtonian part, plus terms due 
to the modification. These in turn can be split into the 
contributions from shells internal and exterior to r. If 
we write the mass overdensity in the ith shell as SMi = 
pip) Si Airrf dri then 



a ■ 



GSM in F z {r s ) - Fi(r c ) 



(A3) 



where 



^(r.) = (l+— ) (e"l r - r *l/ r » - e -< r+r '>/r.) 

= 2 (l + — ) e~ r ' r * smhin/r x ) (A4) 

when r x < r (i.e., from the internal shells) and 

Fi (r x ) = (— - l) e-\ ri - r V r * - (l + — ) e-( r+r *)/ r * 



= -2 



r/r x cosh(r/r x ) - sinh(r/r a ) - Ti / T 
r/r x 



(A5) 



if r x > r (the external shells). Given these expressions, 
we can (and do) integrate to get the total force at a radius 
r. 



To get a qualitative feel for the effect of a ^ 0, it is 
instructive to study the net force associated with a tophat 
perturbation, for which Si is the same for all shells within 
the perturbation. (Note that all the analysis in the main 
text, and the figures in this Appendix, were made by 
solving for the exact evolution, i.e., not assuming that 
the initial tophat remains a tophat.) In this case the 
modification to the Newtonian force —GM/r 2 is given 
by 

F a (r) = a^-[e- r / r °(l + r/r s )T(R/r s ) (A6) 
-e- r / r «(l + r/r c )T(R/r c ) 
1-e-r/r* (i + r /r s )T(R/r s ) . 



GSM 



In the present context, we are also interested in the profile 
within the perturbation. A similar analysis of $(r) when 
r < R shows that the modification to the force is 



GSM{< r) 



- R / r °{l + R/r s )T{r/r s ) 



- R/r ^l + R/r c )T(r/r c ) 



(A7) 



-a- 



GSM{< r) 



l- e - R / r >(l + R/r s )fir/r 3 ) 



The term in square brackets is always positive, so the 
correction is proportional to the product of a and SM . 
Hence, when both are positive, then the resulting profile 
becomes steeper than when a — 0: one might generically 
expect halos to be more concentrated. However, for over- 
dense perturbations when a < 0, the net force towards 
the center is smaller than when a = 0, so we expect the 
profile to evolve away from a tophat in the sense of be- 
coming more concentrated, not at the center, but at the 
outer boundary. 

We have looked for this effect in the simulations 
we study in the main text. Unfortunately, because 
the particle mass in these simulations is rather large 
(the runs were designed to study larger scale struc- 
tures), we can only estimate halo density profiles for 
the rarest most massive objects. In the a — (—1,0,0.5) 
runs, the most massive object had log lo (M//i _1 M ) = 
(14.85,15.15,15.58) and virial radii, ii vir //i" 1 Mpc = 
(1.31,1.66,2.36). Figure [5] shows the profiles of these 
three objects. In each case, we show the expected en- 
closed density profile for an object of this mass if a = 0. 
In all cases, the measured profiles lie below this line on 
scales smaller than i? v i r , suggesting that our simulations 
are not able to properly resolve small scale structures. 
However, there is a hint that the a = — 1 object has 
slightly more than the expected amount of mass just be- 
yond the virial radius, whereas the opposite is true for the 
a = 0.5 object. It will be interesting if these trends per- 
sist in simulations with better mass and force resolution, 
because it is conceivable that effects of this magnitude 
will soon be measured by weak lensing surveys. 
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FIG. 10: One, two, and three a confidence intervals in the 
(a, r s ) plane when the initial power-spectra are the same 
(blue) and when the z = linear theory power spectra are 
the same (black). Constraints come from requiring the mod- 
els to agree with the cluster abundances measured by [4(J , and 
assuming that the scaling relation between cluster observable 
and mass is the same as when a — 0. 

Finally, we note that there is another sense in which 
the modified gravity model is more complicated. As the 
object pulls itself together from the expansion, it will 
eventually reach turnaround. In standard gravity, all its 
shells reach turnaround at the same time, so the total 
energy at turnaround is potential, and energy conser- 
vation means that this equals its initial energy. When 
a ^ 0, however, the different shells turnaround at differ- 
ent times, so the total energy when the outermost shell 
turns around is no longer all potential. 



r s , and an amplitude a. Previous work has sought to con- 
strain the values of these parameters by using the shape 
of the galaxy power spectrum 0,H3|. Although, on large 
scales, the amplitude of the power-spectrum depends on 
a, potentially providing a powerful constraint, this is off- 
set by the fact that galaxies are biased tracers of the dark 
matter field, and this bias is unknown. Even allowing 
for the simplest, linear bias relation removes most of the 
sensitivity of this method to the information contained 
in the amplitude of the late-time P(k). Our succesful 
model for halo abundances, however, is sensitive to pre- 
cisely this difference, so, we can use current high mass 
cluster counts to constrain the allowed range of (a,r s ). 

Figure ITD1 shows the result of comparing our predicted 
mass functions with the measured counts of [40| when we 
normalize to the same (CMB) initial conditions for all a 
(blue contours), and when normalized to have the same 
linear theory P(k) at z = (black). The inner most 
contours represent one, two, and three a confidence lev- 
els. The constraints for the CMB normalized curves are 
slightly weaker than the others. This can be traced to 
Figure 6 in |23j , which shows that the solid line is closer 
to unity than the dotted line in the mass range of clus- 
ters, 10 14 -10 15 M Q //i. This is because the mass function 
is determined by the shape of the barrier and the shape 
of the initial power spectrum, so that in the case of CMB 
normalization only one of these components changes (the 
barrier), whereas for cluster normalization both change, 
leading to a larger change in the mass function - and 
hence to stronger constraints. Note in particular, that 
the constraints that we derive here are roughly compara- 
ble to those derived by [l^, [22J. However, if halo density 
profiles are indeed sensitive to a, as our analysis suggests, 
then this will affect the scaling relations that are usually 
assumed to convert X-ray observables into masses. This, 
in turn, will alter the derived constraints. We look for- 
ward to the day when better simulations are available, 
so the question of how X-ray cluster scaling relations de- 
pend on a and r s has been settled. 



APPENDIX B: CONSTRAINING a AND r s 



The main text studied a particular modification to 
large-scale gravity, which has two free parameters: a scale 



